      Function Kernel(distance)
      include "param.fi"
      include "common.fi"
      Real distance, q, SIGMA
      q = distance / HSUP
      If (KernelType .EQ. WC2) Then
        If (DIMM .EQ. DIMM3D) Then
          SIGMA = 21.0 / 16.0 / PI
          Kernel = (1.0-q)**4.0*(4.0*q+1.0) * SIGMA / H / H / H
        Else
          SIGMA = 7.0 / 4.0 / PI
          Kernel = (1.0-q)**4.0*(4.0*q+1.0) * SIGMA / H / H
        Endif
      Endif
      End Function
C==================================================================
      Subroutine KernelGradient(XA, YA, ZA, XX, YY, ZZ)
      include "param.fi"
      include "common.fi"
      Real distance, q, SIGMA, XA, YA, ZA, XX, YY, ZZ
      distance = SQRT(XA*XA+YA*YA+ZA*ZA)
      q = distance / HSUP
      If (KernelType .EQ. WC2) Then
        If (DIMM .EQ. DIMM3D) Then
          SIGMA = 21.0 / 32.0 / PI
          SIGMA = SIGMA * 20.0 * q * (q-1.0)**3.0
          SIGMA = SIGMA / H / H / H / H / distance
        Else
          SIGMA = 7.0 / 8.0 / PI
          SIGMA = SIGMA * 20.0 * q * (q-1.0)**3.0
          SIGMA = SIGMA / H / H / H / distance
        Endif
        XX = SIGMA * XA
        YY = SIGMA * YA
        ZZ = SIGMA * ZA
      Endif
      End Subroutine
C==================================================================
      Function KernelLaplace(distance)
      include "param.fi"
      include "common.fi"
      Real distance, q, SIGMA
      q = distance / HSUP
      If (KernelType .EQ. WC2) Then
        If (DIMM .EQ. DIMM3D) Then
          SIGMA = 21.0 / 16.0 / PI * 2.0 / 3.0;
          KernelLaplace =
     &    (1.0-q)**2.0*(25.0*q-10.0) * SIGMA / H / H / H / H / H
        Else
          SIGMA = 7.0 / 4.0 / PI * 2.0 / 3.0;
          KernelLaplace =
     &    (1.0-q)**2.0*(25.0*q-10.0) * SIGMA / H / H / H / H
        Endif
      Endif
      End Function
